To conduct a spatio-temporal analysis of tree cover change in a certain country between 2018 and 2023, we obtained geo-spatial data from the Living Atlas. This dataset provides annual land cover information, including tree cover, in a tile-based GeoTiff format. Living Atlas allows users to inspect various types of land cover through satellite images for the time span between 2018 and 2023, including Trees, Water, Crops, Snow/Ice, Built Area, Clouds, Flooded Vegetation, Bare Ground and Rangeland. More details about the data set can be found here.

Since GeoTIFF data of the world is obtained in tiles and is very memory consuming, I decided to start my research based on certain countries of interest, which in our case was Austria. Wile inspecting the data and trying to construct informative visualizations, I ran into some problems with combining rasters that should cover the whole land of Austria. In order to fully resemble Austrian land, I had to download three tiles, but while trying to combine them into one big raster or processing them individually, a lot of data ended up missing, which shouldn’t be the case, since the visualization on the website clearly shows data points in all Austrian areas. We are going to talk about this in second half of this report where I show problems I had with creation of maps for Austria.

To fully test my ideas and see if a comprehensive analysis of GeoTiff data was possible, I decided to start this work with analysis of a country which is fully captured in only one tile - Serbia.

Our first goal is to inspect tree gain and loss in country of Serbia taking years 2018 and 2023 into comparison. Reason we are taking 2018 into account is because Land Cover dataset description pointed out that 2017 (the first year of measurements) shouldn’t be taken into account due to its resemblance of fewer image sets. In order to fully observe whole Serbian land cover property we need to download only one tile which covers the whole aria. We are going to download these tiles for years 2018 and 2023 which we wish to compare:

urls <- c(
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2018/34T_20180101-20190101.tif",
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2018/34T_20180101-20190101.tif"
)

Having the needed download links we are going to programmatically download the tow files and save them into a dedicated folder inside of our working directory:

for (url in urls) {
  download.file(
    url = url,
    destfile = basename("SRB/" + url),
    mode = "wb",
    method = "wget"
  )
}

We are than proceeding with loading necessery libraries which will help us work with country data and GeoTiff images and their corresponding rasters. For those purposes we are going to use library terra (Terra package) for spatial data analysis, sf (Sf package) to work with simple features and conversions of spatial vector data and its projection, exactextractr (Exactextractr package) to accurately summarize raster values over polygonal areas fr zonal statistics and rgeoboundaries (Rgeoboundaries package) which provides access to the geoBoundaries international boundary database here.

For visualizations we are going to use ggplot2, plotly and grid packages.

libs <- c(
  "tidyr", "dplyr", "terra",
  "sf", "exactextractr",
  "rgeoboundaries", "ggplot2",
  "plotly", "grid", "ggtern"
)

installed_libraries <- libs %in% rownames(
  installed.packages()
)

if (any(installed_libraries == F)) {
  install.packages(libs[!installed_libraries])
}

invisible(
  lapply(
    libs, library,
    character.only = T
  )
)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## terra 1.7.78
## 
## Attaching package: 'terra'
## The following object is masked from 'package:tidyr':
## 
##     extract
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
## 
## Attaching package: 'grid'
## The following object is masked from 'package:terra':
## 
##     depth
## Registered S3 methods overwritten by 'ggtern':
##   method           from   
##   grid.draw.ggplot ggplot2
##   plot.ggplot      ggplot2
##   print.ggplot     ggplot2
## --
## Remember to cite, run citation(package = 'ggtern') for further info.
## --
## 
## Attaching package: 'ggtern'
## The following objects are masked from 'package:ggplot2':
## 
##     aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob,
##     ggsave, layer_data, theme_bw, theme_classic, theme_dark,
##     theme_gray, theme_light, theme_linedraw, theme_minimal, theme_void

Loading Country Boundaries of Serbia

For the purposes of this task we are going to load sf objects of both Serbia and Kosovo and bind them. The idea here was to see if this would work as expected and if we would get the right object once both objects are combined. We are going to load administrative bounds of level 1:

serbia_sf <- rgeoboundaries::gb_adm1(
  "SRB"
)
kosovo_sf <- rgeoboundaries::gb_adm1(
  "XKX"
)

We are going to use terra::plot function to see if the to objects were merged correctly:

serbia_sf <- rbind(serbia_sf, kosovo_sf)
terra::plot(serbia_sf)

Loading Raster Files

Next we are going to load the two TIFF files for years 2018 and 2023 into our working environment and convert them into rasters:

raster_files_serbia <- list.files(
  path = "SRB",
  pattern = "tif",
  full.names = T
)

raster_srb_2018 <- terra::rast(raster_files_serbia[1]) 
raster_srb_2023 <- terra::rast(raster_files_serbia[2])

We have to check if our shape file object uses the same Coordinate Reference System (CRS) as the rasters we just loaded. If the two CRS should be different, we are going to transform the Serbia’s shape file to match CRS of the raster objects, so that the two overlap exactly.

crs_raster <- terra::crs(raster_srb_2018)
crs_shapefile <- sf::st_crs(serbia_sf)


crs_raster
## [1] "PROJCRS[\"WGS 84 / UTM zone 34N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 34N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",21,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 18°E and 24°E, northern hemisphere between equator and 84°N, onshore and offshore. Albania. Belarus. Bosnia and Herzegovina. Bulgaria. Central African Republic. Chad. Croatia. Democratic Republic of the Congo (Zaire). Estonia. Finland. Greece. Hungary. Italy. Kosovo. Latvia. Libya. Lithuania. Montenegro. North Macedonia. Norway, including Svalbard and Bjornoys. Poland. Romania. Russian Federation. Serbia. Slovakia. Sudan. Sweden. Ukraine.\"],\n        BBOX[0,18,84,24]],\n    ID[\"EPSG\",32634]]"
crs_shapefile
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
if (crs_raster != crs_shapefile) {
  serbia_sf <- sf::st_transform(serbia_sf, crs = crs_raster)
}

We are going to plot these two to make sure the coordinate systems are the same and that the country’s boundaries can be plotted correctly over the raster object.

terra::plot(raster_srb_2018, main = "Raster and Shapefile Overlay")

if (inherits(serbia_sf, "sf")) {
  serbia_sf_sp <- sf::as_Spatial(serbia_sf)
} else {
  serbia_sf_sp <- serbia_sf
}
terra::plot(serbia_sf, add = TRUE, border = "black", col = adjustcolor("red", alpha.f = 0.5))

Here is also the comparison of the two rasters for the years 2018 and 2023:

terra::plot(raster_srb_2018)

terra::plot(raster_srb_2023)

Extracting meaningful data from rasters

In following we are going to take rasters for the two years (2018 and 2023) and try to extract meaningful data from them regarding tree cover change that happened in span of those 5 years. Besides that we also want to capture built areas of the country and changes in area covered by them in that time as well.

First we are going to make data frame out of the color table of our raster to inspect with which colors we are working since the colors represent certain values in our raster - meaning certain color at one point stands for either built area, a tree or something else that the data set captures. Here we want to concentrate us on colors representing trees and built areas and values associated with them so that we can use these for further visualization purposes. To get the idea about the colors we are going to use terra::coltab function:

raster_color_table <- do.call(
    data.frame,
    terra::coltab(raster_srb_2023)
)
head(raster_color_table)

Taking a look at the first few rows of the color palette we can see that each color value is separated in its RGB features and captured over columns 2,3,4 featuring also the alpha column which is not meaningful for us. First column named value represents the value through which we can access that color in the data set, which also means that through that value of the raster point we can access only points that are meaningful for us.

Since we work with hexadecimal codes for colors in R we are going to convert these into their hexadecimal representatives using ggtern::rgb2hex function.

hex_code <- ggtern::rgb2hex(
    r = raster_color_table[,2],
    g = raster_color_table[,3],
    b = raster_color_table[,4]
)
hex_code
##   [1] "#000000" "#419bdf" "#397d49" "#000000" "#7a87c6" "#e49635" "#000000"
##   [8] "#c4281b" "#a59b8f" "#a8ebff" "#616161" "#e3e2c3" "#000000" "#000000"
##  [15] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [22] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [29] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [36] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [43] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [50] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [57] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [64] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [71] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [78] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [85] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [92] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
##  [99] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [106] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [113] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [120] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [127] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [134] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [141] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [148] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [155] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [162] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [169] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [176] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [183] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [190] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [197] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [204] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [211] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [218] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [225] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [232] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [239] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [246] "#000000" "#000000" "#000000" "#000000" "#000000" "#000000" "#000000"
## [253] "#000000" "#000000" "#000000" "#000000"

Taking a look at hexadecimal codes we see that only first 12 entries hold actual colors. Everything else has code #000000 which is black and represents a missing value. So the colors we get are:

c("#000000", "#419bdf", "#397d49", "#000000", "#7a87c6", "#e49635", "#000000", "#c4281b", "#a59b8f", "#a8ebff" ,"#616161" ,"#e3e2c3")
##  [1] "#000000" "#419bdf" "#397d49" "#000000" "#7a87c6" "#e49635" "#000000"
##  [8] "#c4281b" "#a59b8f" "#a8ebff" "#616161" "#e3e2c3"

Since we want to extract only trees and built areas which are shown to be green and red on Living Atlas, we want to consider only values with colors coded through #397d49 and #c4281b. We are going to convert this hexadecimal value back to its RGB values using function col2rgb, so that we can find indexes of colors which we will be extracting, in our case green for trees and red for built areas:

built_target_color <- "#c4281b"
built_target_rgb <- col2rgb(built_target_color)

tree_target_color <- "#397d49"
tree_target_rgb <- col2rgb(tree_target_color)

We are now going to write a function which is going to help us find the color matches in the color table and extract the values associated with them through which we can later on access the raster point:

compare_rgb <- function(r, g, b, target_rgb) {
  return(r == target_rgb[1] & g == target_rgb[2] & b == target_rgb[3])
}

We will now add two new columns to the raster_color_table: match_built and match_tree. For each row this function will assign TRUE if the color was successfully found and FALSE otherwise.

raster_color_table$match_built <- with(raster_color_table, 
                                       compare_rgb(red, green, blue, built_target_rgb))

raster_color_table$match_tree <- with(raster_color_table, 
                                      compare_rgb(red, green, blue, tree_target_rgb))

Now we can extract the two values and use them in further work for accessing tree and built points:

built_target_index <- raster_color_table$value[raster_color_table$match_built]
tree_target_index <- raster_color_table$value[raster_color_table$match_tree]

We will construct a list of two data frames, first for the year 2018 and the second one fro 2023. Using exactextractr::exact_extract function, we summarize the area covered with the trees and the area covered with cities or infrastructure (built area) for both years. Since coverage_area returns m2 we are going to divide returned values through 1e6 to get km2. To extract tree cover points we are going to consider our previously found tree_target_index and to get information about built areas built_target_index

rasters_srb_both_years <- list(raster_srb_2018, raster_srb_2023)

dd_serbia <- list()
for (i in seq_along(rasters_srb_both_years)) {
  rasters <- rasters_srb_both_years[[i]]
  lc <- exactextractr::exact_extract(
    rasters,
    serbia_sf,
    function(df) {
      df |>
        dplyr::group_by(
          shapeName
        ) |>
        dplyr::summarize(
          tree_area_km2 = sum(
            coverage_area[value == tree_target_index] / 1e6
          ),
          built_area_km2 = sum(
            coverage_area[value == built_target_index] / 1e6
          )
        )
    },
    summarize_df = TRUE,
    coverage_area = TRUE,
    include_cols = "shapeName"
  )
  dd_serbia[[i]] <- lc
}
## Cannot preload entire working area of 1603087171 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.
##   |                                                                              |                                                                      |   0%  |                                                                              |==                                                                    |   3%  |                                                                              |====                                                                  |   6%  |                                                                              |=======                                                               |   9%  |                                                                              |=========                                                             |  12%  |                                                                              |===========                                                           |  16%  |                                                                              |=============                                                         |  19%  |                                                                              |===============                                                       |  22%  |                                                                              |==================                                                    |  25%  |                                                                              |====================                                                  |  28%  |                                                                              |======================                                                |  31%  |                                                                              |========================                                              |  34%  |                                                                              |==========================                                            |  38%  |                                                                              |============================                                          |  41%  |                                                                              |===============================                                       |  44%  |                                                                              |=================================                                     |  47%  |                                                                              |===================================                                   |  50%  |                                                                              |=====================================                                 |  53%  |                                                                              |=======================================                               |  56%  |                                                                              |==========================================                            |  59%  |                                                                              |============================================                          |  62%  |                                                                              |==============================================                        |  66%  |                                                                              |================================================                      |  69%  |                                                                              |==================================================                    |  72%  |                                                                              |====================================================                  |  75%  |                                                                              |=======================================================               |  78%  |                                                                              |=========================================================             |  81%  |                                                                              |===========================================================           |  84%  |                                                                              |=============================================================         |  88%  |                                                                              |===============================================================       |  91%  |                                                                              |==================================================================    |  94%  |                                                                              |====================================================================  |  97%  |                                                                              |======================================================================| 100%
## Cannot preload entire working area of 1603087171 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.
##   |                                                                              |                                                                      |   0%  |                                                                              |==                                                                    |   3%  |                                                                              |====                                                                  |   6%  |                                                                              |=======                                                               |   9%  |                                                                              |=========                                                             |  12%  |                                                                              |===========                                                           |  16%  |                                                                              |=============                                                         |  19%  |                                                                              |===============                                                       |  22%  |                                                                              |==================                                                    |  25%  |                                                                              |====================                                                  |  28%  |                                                                              |======================                                                |  31%  |                                                                              |========================                                              |  34%  |                                                                              |==========================                                            |  38%  |                                                                              |============================                                          |  41%  |                                                                              |===============================                                       |  44%  |                                                                              |=================================                                     |  47%  |                                                                              |===================================                                   |  50%  |                                                                              |=====================================                                 |  53%  |                                                                              |=======================================                               |  56%  |                                                                              |==========================================                            |  59%  |                                                                              |============================================                          |  62%  |                                                                              |==============================================                        |  66%  |                                                                              |================================================                      |  69%  |                                                                              |==================================================                    |  72%  |                                                                              |====================================================                  |  75%  |                                                                              |=======================================================               |  78%  |                                                                              |=========================================================             |  81%  |                                                                              |===========================================================           |  84%  |                                                                              |=============================================================         |  88%  |                                                                              |===============================================================       |  91%  |                                                                              |==================================================================    |  94%  |                                                                              |====================================================================  |  97%  |                                                                              |======================================================================| 100%

Now we are going to bind the columns and get a full data set containing important information for both years:

output <- do.call(cbind, dd_serbia)
serbia_df <- as.data.frame(output)
head(serbia_df)

Renaming columns for clarity:

names(serbia_df) <- c("shapeName", "tree_area_km2_2018", "built_area_km2_2018", "shapeName2",
                      "tree_area_km2_2023", "built_area_km2_2023")

Since we have shapeName column twice, we are going to extract only meaningful columns:

serbia_df <- serbia_df |> dplyr::select(1:3, 5:6)

Last thing left to do is to add two new columns to our data set tree_cover_change and bult_area_change which capture the percentage change in tree cover and built areas in 2018 and 2023.

serbia_df <- serbia_df |>
  dplyr::mutate(
    tree_cover_change = (
      tree_area_km2_2023 - tree_area_km2_2018
      ) / tree_area_km2_2018 * 100
  )

serbia_df <- serbia_df |>
  dplyr::mutate(
    bult_area_change = (
      built_area_km2_2023 - built_area_km2_2018
      ) / built_area_km2_2018 * 100
  )

Lastly we have to bind columns of country’s shape file and created data set, so that this can be used for plotting, because of geometry column in serbia_sf which holds information about polygon to which we are referring as administrative bound of level 1:

serbia_changes <- cbind(serbia_sf, serbia_df)

Visualization

Map of Tree Cover Change in Serbia comparing years 2018 and 2023

map <- ggplot() +
  geom_sf(
    data = serbia_changes,
      aes(
      fill = tree_cover_change),
    color = "white",
    size = .15
  ) +
  scale_fill_gradient2(
    name = "Tree cover change in Serbia comparing 2018 and 2023",
    midpoint = 0,
    mid = "#f7de7c",
    high = "#006f00",
    low = "#9e319d"
  ) +
  guides(
    fill = guide_colorbar(
      direction = "horizontal",
      barheight = unit(1.5, "mm"),
      barwidth = unit(20, "mm"),
      title.position = "top",
      label.position = "bottom",
      title.hjust = .5,
      label.hjust = .5,
      nrow = 1,
      byrow = T
    )
  ) +
  theme_void() +
  theme(
    legend.position = "top",
    plot.margin = unit(
      c(
        t = 0, b = 0,
        r = 0, l = 0
      ), "lines"
    )
  )
map

Map of Built Area Change in Serbia comparing years 2018 and 2023

built_map_change <- ggplot() +
  geom_sf(
    data = serbia_changes,
    aes(
      fill = bult_area_change
    ),    color = "white",
    size = .15
  ) +
  scale_fill_gradient2(
    name = "Built AreaChange in Serbia comparing 2018 and 2023",
    midpoint = 0,
    mid = "yellow",
    high = "red",
    low = "green"
  ) +
  guides(
    fill = guide_colorbar(
      direction = "horizontal",
      barheight = unit(1.5, "mm"),
      barwidth = unit(20, "mm"),
      title.position = "top",
      label.position = "bottom",
      title.hjust = .5,
      label.hjust = .5,
      nrow = 1,
      byrow = T
    )
  ) +
  theme_void() +
  theme(
    legend.position = "top",
    plot.margin = unit(
      c(
        t = 0, b = 0,
        r = 0, l = 0
      ), "lines"
    )
  )
built_map_change

Adding Built Areas To The map of Tree Cover Change

The idea of this section is to plot points of country’s built areas over the the previously constructed map for tree cover change. To do this we will use the data for year 2023. This section is meant for visualization purposes only and further inspections of the relation between tree cover change and gain in built areas are planed. It would be meaningful to perform statistic analysis on this variables and see how they correlate, but in this submission I will be concentrating on inspecting my possibilities with raster data and libraries used for working with them.

Extracting Built Area Values

Function terra::ifel lets us traverse the the raster_srb_2023 based on certain condition and keeping values that we are looking for while maksing out everything else found as NA's. We will now make a use of our variable built_target_index, which represents raster points of built areas and to those points we will assign value 1 and mask everything else out.

built_up_raster <- terra::ifel(
    raster_srb_2023 == built_target_index,
    1, 
    NA 
)
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          

Taking a look at the current plot to see if the extraction was successful:

terra::plot(built_up_raster)

Since built_up_raster also contains points for built areas in neighbouring countries included in initial raster_srb_2023, we will use terra::crop function to consider only the points obtained inside of Serbia’s boundaries given by shape file named serbia_sf. Since terra::crop allows only cropping only rectangular areas we are left with few more points that should be masked out and not considered, so we are going to use function terra::mask to do so:

cropped_raster <- crop(built_up_raster, vect(serbia_sf), snap = "in")
## |---------|---------|---------|---------|=========================================                                          
masked_raster <- mask(cropped_raster, vect(serbia_sf))
## |---------|---------|---------|---------|=========================================                                          

Plotting masked_raster to see if it is sufficient:

terra::plot(masked_raster)

We will now construct a data frame from our masked_raster, so that we can really plot these points over our tree cover change map. We are also going to use aggregate function with factor 5, so that memory consumption stays under control since we are working with huge rasters:

built_up_df <- as.data.frame(aggregate(masked_raster, fact = 5), xy = TRUE, na.rm = TRUE)
## |---------|---------|---------|---------|=========================================                                          
head(built_up_df)

The third column of our new data frame has always the value 1, which means that each coordinate in the raster represents a point that stands for a built area. We are going to convert this column to a factor, so that no problems with plotting occure:

names(built_up_df)[3] <- "layer"
built_up_df <- built_up_df %>%
  mutate(layer = as.factor(layer))

Plotting the map of built areas

Let us see how the map of built areas is looking like so far:

built_map <- ggplot() +
  geom_sf(data = serbia_sf, fill = NA, color = "black") +
  geom_tile(data = built_up_df, aes(x = x, y = y, fill = layer), show.legend = FALSE) +
  scale_fill_manual(values = built_target_color) +
  theme_minimal() +
  labs(title = "Built-Up Areas in Serbia (2023)")
built_map

Since we got what we initially wanted, we can make a combined map of tree cover change overlayed by map of built areas:

built_up_color <- "red"

combined_map <- ggplot() +
  geom_sf(data = serbia_changes, aes(fill = tree_cover_change), color = "white", size = 0.15) +
  scale_fill_gradient2(
    name = "Tree cover change in Serbia comparing 2018 and 2023",
    midpoint = 0,
    mid = "#f7de7c",
    high = "#006f00",
    low = "#9e319d",
    guide = guide_colorbar(
      direction = "horizontal",
      barheight = unit(1.5, "mm"),
      barwidth = unit(20, "mm"),
      title.position = "top",
      label.position = "bottom",
      title.hjust = 0.5,
      label.hjust = 0.5,
      nrow = 1,
      byrow = TRUE
    )
  ) +
  # Overlay built-up areas with a fixed color
  geom_tile(data = built_up_df, aes(x = x, y = y), fill = built_up_color, show.legend = FALSE) +
  theme_void() +
  theme(
    legend.position = "top",
    plot.margin = unit(c(t = 0, b = 0, r = 0, l = 0), "lines"),
    plot.title = element_text(hjust = 0.5)  # Center the title
  ) +
  labs(title = "Built-Up Areas in Serbia (2023)")

combined_map

Adding City To The Map

To further enhance our map we are going to add city names to it to represent the biggest cities in the country (which are also the biggest built areas). To do so we are going to download data set Geonames Cities with population > 5000 which contains all relevant information of cities including population of the city, country code and geographical location of the city which we need for plotting:

filename <- "geonames-population.csv"
# download
get_geonames_data <- function(){
    table_link <- "https://documentation-resources.opendatasoft.com/api/explore/v2.1/catalog/datasets/doc-geonames-cities-5000/exports/csv?lang=en&timezone=Europe%2FBerlin&use_labels=true&delimiter=%3B"
    res <- httr::GET(
        table_link,
        httr::write_disk(
            filename
        ),
        httr::progress()
    )
}

get_geonames_data()

Let’s now read the data set we just downloaded:

load_geonames_data <- function(){
    places_df <- read.csv(
        filename,
        sep = ";"
    )
    return(places_df)
}

places_df <- load_geonames_data()

head(places_df)

Columns we are interested in are 2 (name), 7 (country_code), 13 (population) and 18 (location). We are now going to store these in a separate data frame named places_modified_df. Since column location has values in form latitude, longitude we are going to split these values into two separate columns named lat and long:

places_modified_df <- places_df[, c(2, 7, 13, 18)]
names(places_modified_df) <- c(
    "name", "country_code", "pop", "coords")

places_modified_df[c("lat", "long")] <-
    stringr::str_split_fixed(
        places_modified_df$coords, ",", 2
    )

Now we have to make a shape file which we are going to be able to plot. This is a for of data frame where we will extract only cities with country_code “RS” or “XK” to fit to our map and based on population variable we are going to take the 5 biggest cities. To make it a real sf object we also need coordinates which are constructed from our variables long and lat. Lastly we have to project this in coordinate system used for Europe which is: "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" :

crs_lambert <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
places_clean_sf <- places_modified_df |>
    dplyr::filter(country_code == "RS" | country_code == "XK") |>
    dplyr::slice_max(
        pop,
        n = 5
    ) |>
    dplyr::select(
        -coords,
        -country_code,
        -pop
    ) |>
    sf::st_as_sf(
        coords = c(
            "long",
            "lat"
        ),
        crs = 4326
    ) |>
    sf::st_transform(crs = crs_lambert)

Our complete map looks like this now:

cities_added <- combined_map + 
    geom_sf_text(
    data = places_clean_sf,
    aes(label = name),
    size = 3,      
    color = "black" 
  )
cities_added

Austria

I initially tried to perform these visualizations on rasters for Austria, but as I’ve mentioned at the beginning of the document I ran into some problems when merging 3 raster object that should represent entire country of Austria.

First lets download all GeoTIF files for Austria for years 2018 and 2023:

urls <- c(
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2023/33T_20230101-20240101.tif",
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2018/33T_20180101-20190101.tif",
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2023/33U_20230101-20240101.tif",
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2018/33U_20180101-20190101.tif",
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2023/32T_20230101-20240101.tif",
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2018/32T_20180101-20190101.tif"
)
for (url in urls) {
  download.file(
    url = url,
    destfile = basename("AT/" + url),
    mode = "wb",
    method = "wget"
  )
}

Administrative Boundary Data

The analysis was focused on Austria, for which administrative boundary data were obtained using the rgeoboundaries package. Specifically, second-level administrative boundaries (e.g., districts or municipalities) were retrieved:

austria_sf <- rgeoboundaries::gb_adm1(
  "AUT"
)

Plot Austria:

terra::plot(austria_sf)

Raster Data Processing

I than rearranged files for the two years in two separate folders inside of AT directory, one folder 2018 and one 2023 so the access is easier. Raster files corresponding to the tree cover data for 2018 and 2023 were identified and loaded using the terra package. The files were selected based on the .tif file extension, ensuring that all relevant data were included:

raster_files_2018 <- list.files(
  path = "AT/2018",
  pattern = "20190101.tif",
  full.names = T
)
raster_files_2023 <- list.files(
  path = "AT/2023",
  pattern = "20240101.tif",
  full.names = T
)

My idea than was to merge all GeoTIF files for one year into one single raster, so that I get a full image at once, that can be accomplished by first converting all of the GeoTIF files into rasters, than using sprc to create collection of the rasters from the list of rasters which than can be combined to a single raster using terra::mosaic function:

rasters_2018 <- lapply(raster_files_2018, terra::rast)
rasters_combined_2018 <- terra::mosaic(sprc(rasters_2018), fun = mean)
## |---------|---------|---------|---------|=========================================                                          

As for Serbia I tried the same approach with converting the shape file and the raster to the same CRS:

crs_raster_aut <- terra::crs(rasters_combined_2018)
crs_shapefile_aut <- sf::st_crs(austria_sf)


crs_raster_aut
## [1] "PROJCRS[\"WGS 84 / UTM zone 32N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 32N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n        BBOX[0,6,84,12]],\n    ID[\"EPSG\",32632]]"
crs_shapefile_aut
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
if (crs_raster_aut != crs_shapefile_aut) {
  austria_sf <- sf::st_transform(austria_sf, crs = crs_raster_aut)
}

But than as soon as I plotted my raster only a small part of Austria was shown on the map, seemingly not at the right place at all:

terra::plot(rasters_combined_2018, main = "Raster and Shapefile Overlay Austria")

if (inherits(austria_sf, "sf")) {
  austria_sf_sp <- sf::as_Spatial(austria_sf)
} else {
  austria_sf_sp <- austria_sf
}
terra::plot(austria_sf, add = TRUE, border = "black", col = adjustcolor("red", alpha.f = 0.5))

I was wandering if something went wrong while merging rasters using mosaic function, so I tried to plot at least retrospective parts of Austrian maps on each of three rasters:

for (rast in rasters_2018) {
  terra::plot(rast, main = "Tiles of austria raster and overlay of Austrian borders")

if (inherits(austria_sf, "sf")) {
  austria_sf_sp <- sf::as_Spatial(austria_sf)
} else {
  austria_sf_sp <- austria_sf
}
terra::plot(austria_sf, add = TRUE, border = "black", col = adjustcolor("red", alpha.f = 0.5))
}

As it can be seen it appears that the rasters weren’t combined successfully into one. North Sea isn’t shown on the combined raster as well as Sicily which are seen on the separate tiles. My next thought was that the resulting raster maybe has to be in rectangular shape in order to be drawn properly, although I didn’t see anything about this in the terra documentation. Now we will download one more GeoTIF file, so that we get a complete rectangle around Austria:

missing_tile <- "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2018/32U_20180101-20190101.tif"
destfile_path <- paste0("AT/2018/MISSING_TILE/", basename(missing_tile))

Download the tile:

download.file(
  url = missing_tile,
  destfile = destfile_path,
  mode = "wb",
  method = "wget"
)

Now I will add this raster to rasters_2018 and try merging them into one single raster again:

missing_raster <- terra::rast(destfile_path)
rasters_2018 <- append(rasters_2018, missing_raster)

Now lets mosaic all four rasters together:

rasters_combined_2018 <- terra::mosaic(sprc(rasters_2018), fun = mean)
## |---------|---------|---------|---------|=========================================                                          

Lets plot it to see how it looks like now:

terra::plot(rasters_combined_2018, main = "Rasters capturing Austria")

Again the same plot is shown as result of mosaicing the four tiles. Since I assumed all rasters have the same CRS, I decided to check this and allign all of the rasters to CRS of the first one:

lapply(rasters_2018, terra::crs)
## [[1]]
## [1] "PROJCRS[\"WGS 84 / UTM zone 32N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 32N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n        BBOX[0,6,84,12]],\n    ID[\"EPSG\",32632]]"
## 
## [[2]]
## [1] "PROJCRS[\"WGS 84 / UTM zone 33N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 33N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",15,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State.\"],\n        BBOX[0,12,84,18]],\n    ID[\"EPSG\",32633]]"
## 
## [[3]]
## [1] "PROJCRS[\"WGS 84 / UTM zone 33N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 33N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",15,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State.\"],\n        BBOX[0,12,84,18]],\n    ID[\"EPSG\",32633]]"
## 
## [[4]]
## [1] "PROJCRS[\"WGS 84 / UTM zone 32N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 32N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.\"],\n        BBOX[0,6,84,12]],\n    ID[\"EPSG\",32632]]"

To make sure all rasters are in same CRS I manually assigned them coordinate reference system associated with Central Europe EPSG:32632:

target_crs <- "EPSG:32632"

rasters_2018 <- lapply(rasters_2018, function(r) {
  if (is.null(terra::crs(r))) {
    terra::crs(r) <- target_crs
  }
  r
})

Than I also checked the resolutions of the rasters and resampled other rasters to resolution of the first:

common_res <- terra::res(rasters_2018[[1]])
rasters_2018 <- lapply(rasters_2018, function(r) {
  if (!all(terra::res(r) == common_res)) {
    r <- terra::resample(r, rasters_2018[[1]], method = "bilinear")
  }
  r
})

I also checked extents of the rasters to make sure they are the same:

common_extent <- terra::ext(rasters_2018[[1]])
rasters_2018 <- lapply(rasters_2018, function(r) {
  if (!identical(terra::ext(r), common_extent)) {
    r <- terra::extend(r, common_extent)
  }
  r
})
## |---------|---------|---------|---------|=================================                                          |---------|---------|---------|---------|=========================================                                          

Again I tried to mosaic the rasters together, but this time aggregating them, so the function finishes faster:

rasters_combined_2018 <- terra::mosaic(terra::sprc(lapply(rasters_2018, function (r) {
  aggregate(r,
            fact = 10)
})), fun = mean)
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          

The end result was again the same and not what I expected to see:

terra::plot(rasters_combined_2018, main = "Merged Raster of Austria 2018")

Using Virtual Layer To Join Rasters

We concluded that each of the tree rasters has diferent CRS and that terra::mosaic isn’t able to combine all tree raster files to one single raster. I researched and found that the standard CRS is EPSG:4326. Since all tree rasters have different CRS, that is the reason why Austrian map can’t be converted to one of them and than be used for all of them repeatedly. in following portion of code, I’ll transform austria_sf to the CRS of the according raster file. Than I will be able to crop the right part of Austrian land from the raster using functions from terra package crop and mask. The raster is than projected to Coordinate reference system EPSG:4326. Than those new rasters will be saved into our directory AUT/2018.

#remove missing raster
rasters_2018[[4]] <- NULL
i <- 1
crs <- "EPSG:4326"
for(raster in rasters_2018){
  country <- austria_sf %>%
        sf::st_transform(crs = terra::crs(raster))
  
 land_cover <- crop(raster, vect(country), snap = "in") %>%
    mask(vect(country)) %>%
    aggregate(fact = 5, fun = "modal") %>%
    terra::project(crs)

    terra::writeRaster(
        land_cover,
        paste0("AT/2018/",
            i,
            "_austria_2018",
            ".tif"
        ),
        overwrite = TRUE
    )
    i <- i + 1
}
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          

Having our new rasters, we can now load them into our session as a list and use them to make a virtual layer, which is basically the same as combining those rasters using mosaic function:

r_list <- list.files(
    path = "AT/2018",
    pattern = "_austria_2018",
    full.names = T
)
land_cover_vrt_2018 <- terra::vrt(
    r_list,
    "austria_land_cover_vrt_2018.vrt",
    overwrite = T
)

Plotting the new virtual land cover raster, revealed, that we achieved the desired behaviour and now we have complete Austrian map with all points from the Land Cover data set. The only difference is that the color of the map are not the same as original ones after cropping and masking of only Austrian land was performed:

terra::plot(land_cover_vrt_2018)

We will now do the same for rasters from 2023, first by rasterizing the GeoTIF files:

rasters_2023 <- lapply(raster_files_2023, terra::rast)

Now let us project the the rasters to appropriate CRS and save them:

i <- 1
for(raster in rasters_2023){
  country <- austria_sf %>%
        sf::st_transform(crs = terra::crs(raster))
  land_cover <- crop(raster, vect(country), snap = "in") %>%
    mask(vect(country)) %>%
    aggregate(fact = 5, fun = "modal") %>%
    terra::project(crs)

    terra::writeRaster(
        land_cover,
        paste0("AT/2023/",
            i,
            "_austria_2023",
            ".tif"
        ),
        overwrite = TRUE
    )
    i <- i + 1
}
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          

Constructing virtual layer for year 2023:

r_list_23 <- list.files(
    path = "AT/2023",
    pattern = "_austria_2023",
    full.names = T
)
land_cover_vrt_2023 <- terra::vrt(
    r_list_23,
    "austria_land_cover_vrt_2023.vrt",
    overwrite = T
)

The plot:

terra::plot(land_cover_vrt_2023)

Retrieving Original Colors

I will once again inspect raster colors. Since we are using instances of Land Cover Data Set, the colors should be same, as they were for Serbia in previous part of this submission, but we will check one more time to be sure. For that it is enough to just take a look at one original raster:

raster_color_table_aut <- do.call(
    data.frame,
    terra::coltab(rasters_2018[[1]])
)

head(raster_color_table_aut)

As we supposed the colors are the same. This time we will keep all colors for visualization purposes:

cols <- hex_code[c(2:3, 5:6, 8:12)]

from <- c(1:2, 4:5, 7:11)
to <- t(col2rgb(cols))
land_cover_vrt_2018 <- na.omit(land_cover_vrt_2018)

To change colors in our raster we will substitute them using terra::subst function:

land_cover_at_2018 <- terra::subst(
    land_cover_vrt_2018,
    from = from,
    to = to,
    names = cols
)
## |---------|---------|---------|---------|=========================================                                          

We will do the same for year 2023:

land_cover_at_2023 <- terra::subst(
    land_cover_vrt_2023,
    from = from,
    to = to,
    names = cols
)
## |---------|---------|---------|---------|=========================================                                          

Plotting my previously created rasters land_cover_at_2018 and land_cover_at_2023, I realized that those are multi-band rasters, meaning red, green and blue are divided into separate bands, which would require using terra::plotRGB function and would make extracting single colors more difficult:

terra::plot(land_cover_at_2018)

terra::plot(land_cover_at_2023)

Plotting both next to each other using terra::plotRGB:

terra::plotRGB(land_cover_at_2018, main = "Land Cover in Austria 2018")
legend("topleft", legend = c("Water", "Trees", 
                             "Flooded vegetation", "Crops", "Built Area",
                             "Bare ground",
                             "Snow/Ice",
                             "Clouds", "Rangeland"), pch = 20, xpd=NA, cex = 0.6, bg="white", 
       col=c("#419bdf", "#397d49", "#7a87c6", "#e49635", 
             "#c4281b", "#a59b8f", "#a8ebff" ,"#616161" ,"#e3e2c3"))

terra::plotRGB(land_cover_at_2023, main = "Land Cover in Austria 2023")
legend("topleft", legend = c("Water", "Trees", 
                             "Flooded vegetation", "Crops", "Built Area",
                             "Bare ground",
                             "Snow/Ice",
                             "Clouds", "Rangeland"), pch = 20, xpd=NA, cex = 0.6, bg="white", 
       col=c("#419bdf", "#397d49", "#7a87c6", "#e49635", 
             "#c4281b", "#a59b8f", "#a8ebff" ,"#616161" ,"#e3e2c3"))

Extracting Data from Virtual Rasters

In order to come up with the same maps as we did for Serbia, we will proceed with extracting values from the two virtual rasters, as they are the ones still containing values needed for point extraction:

virtual_rasters <- list(land_cover_vrt_2018, land_cover_vrt_2023)

dd_austria <- list()
for (i in seq_along(virtual_rasters)) {
  rasters <- virtual_rasters[[i]]
  lc <- exactextractr::exact_extract(
    rasters,
    austria_sf %>% sf::st_transform(crs = terra::crs(rasters)),
    function(df) {
      df %>%
        dplyr::group_by(
          shapeName
        ) %>%
        dplyr::summarize(
          tree_area_km2 = sum(
            coverage_area[value == tree_target_index] / 1e6
          ),
          built_area_km2 = sum(
            coverage_area[value == built_target_index] / 1e6
          )
        )
    },
    summarize_df = TRUE,
    coverage_area = TRUE,
    include_cols = "shapeName"
  )
  dd_austria[[i]] <- lc
}
## Cannot preload entire working area of 54919914 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.
##   |                                                                              |                                                                      |   0%  |                                                                              |========                                                              |  11%  |                                                                              |================                                                      |  22%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================                                       |  44%  |                                                                              |=======================================                               |  56%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================                |  78%  |                                                                              |==============================================================        |  89%  |                                                                              |======================================================================| 100%
## Cannot preload entire working area of 51625668 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.
##   |                                                                              |                                                                      |   0%  |                                                                              |========                                                              |  11%  |                                                                              |================                                                      |  22%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================                                       |  44%  |                                                                              |=======================================                               |  56%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================                |  78%  |                                                                              |==============================================================        |  89%  |                                                                              |======================================================================| 100%
dd_austria[[1]]
dd_austria[[2]]

Both resulting data sets include only values for Vienna, which is concerning since there should be data available for each region. I decided to overlay land_cover_vrt_2018 with austria_sf to see if the match correctly, which is true:

terra::plot(land_cover_vrt_2018, main = "Raster and Shapefile Overlay Austria")
terra::plot(austria_sf  %>% sf::st_transform(crs = terra::crs(rasters)),
            add = TRUE, border = "black", col = adjustcolor("red", alpha.f = 0.3))

I than decided to make a simple text extraction to see if the data would still be missing:

single_region <- austria_sf %>% 
  filter(shapeName == "Burgenland") %>%
  sf::st_transform(crs = terra::crs(land_cover_vrt_2018))  

test_extract <- exactextractr::exact_extract(
  land_cover_vrt_2018,
  single_region,
  function(df) {
    print(head(df))  
    df %>%
      dplyr::summarize(
        tree_area_km2 = sum(
          coverage_area[value == tree_target_index] / 1e6,
          na.rm = TRUE
        ),
        built_area_km2 = sum(
          coverage_area[value == built_target_index] / 1e6,
          na.rm = TRUE
        )
      )
  },
  summarize_df = TRUE,
  coverage_area = TRUE,
  include_cols = "shapeName"
)
##   value  shapeName coverage_area
## 1   NaN Burgenland     424.07799
## 2   NaN Burgenland      79.02509
## 3     5 Burgenland     611.78419
## 4     5 Burgenland    2873.86222
## 5   NaN Burgenland    2775.98223
## 6   NaN Burgenland    1441.25257
test_extract

The inspection has shown, that I was forgetting to exclude missing values (NAs) while summarizing the values for each province:

virtual_rasters <- list(land_cover_vrt_2018, land_cover_vrt_2023)

dd_austria <- list()
for (i in seq_along(virtual_rasters)) {
  rasters <- virtual_rasters[[i]]
  lc <- exactextractr::exact_extract(
    rasters,
    austria_sf %>% sf::st_transform(crs = terra::crs(rasters)),
    function(df) {
      df %>%
        dplyr::group_by(
          shapeName
        ) %>%
        dplyr::summarize(
          tree_area_km2 = sum(
            coverage_area[value == tree_target_index] / 1e6,
             na.rm = TRUE

          ),
          built_area_km2 = sum(
            coverage_area[value == built_target_index] / 1e6,
            na.rm = TRUE
          )
        )
    },
    summarize_df = TRUE,
    coverage_area = TRUE,
    include_cols = "shapeName"
  )
  dd_austria[[i]] <- lc
}
## Cannot preload entire working area of 54919914 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.
##   |                                                                              |                                                                      |   0%  |                                                                              |========                                                              |  11%  |                                                                              |================                                                      |  22%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================                                       |  44%  |                                                                              |=======================================                               |  56%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================                |  78%  |                                                                              |==============================================================        |  89%  |                                                                              |======================================================================| 100%
## Cannot preload entire working area of 51625668 cells with max_cells_in_memory = 3e+07. Raster values will be read for each feature individually.
##   |                                                                              |                                                                      |   0%  |                                                                              |========                                                              |  11%  |                                                                              |================                                                      |  22%  |                                                                              |=======================                                               |  33%  |                                                                              |===============================                                       |  44%  |                                                                              |=======================================                               |  56%  |                                                                              |===============================================                       |  67%  |                                                                              |======================================================                |  78%  |                                                                              |==============================================================        |  89%  |                                                                              |======================================================================| 100%

Now we can move forward and combine the data into one single data frame ready for plotting:

output <- do.call(cbind, dd_austria)
austria_df <- as.data.frame(output)
head(austria_df)
names(austria_df) <- c("shapeName", "tree_area_km2_2018", "built_area_km2_2018", "shapeName2",
                      "tree_area_km2_2023", "built_area_km2_2023")

austria_df <- austria_df |> dplyr::select(1:3, 5:6)
austria_df <- austria_df |>
  dplyr::mutate(
    tree_cover_change = (
      tree_area_km2_2023 - tree_area_km2_2018
      ) / tree_area_km2_2018 * 100
  )

austria_df <- austria_df |>
  dplyr::mutate(
    bult_area_change = (
      built_area_km2_2023 - built_area_km2_2018
      ) / built_area_km2_2018 * 100
  )
austria_changes <- cbind(austria_sf, austria_df)

Austria Tree Cover Change Map

tree_cover_change_map_austria <- ggplot() +
  geom_sf(
    data = austria_changes,
      aes(
      fill = tree_cover_change),
    color = "white",
    size = .15
  ) +
  scale_fill_gradient2(
    name = "Tree cover change in Austria comparing 2018 and 2023",
    midpoint = 0,
    mid = "#f7de7c",
    high = "#006f00",
    low = "#9e319d"
  ) +
  guides(
    fill = guide_colorbar(
      direction = "horizontal",
      barheight = unit(1.5, "mm"),
      barwidth = unit(20, "mm"),
      title.position = "top",
      label.position = "bottom",
      title.hjust = .5,
      label.hjust = .5,
      nrow = 1,
      byrow = T
    )
  ) +
  theme_void() +
  theme(
    legend.position = "top",
    plot.margin = unit(
      c(
        t = 0, b = 0,
        r = 0, l = 0
      ), "lines"
    )
  )
tree_cover_change_map_austria

Austria Built Area Cover Change Map

built_map_change_austria <- ggplot() +
  geom_sf(
    data = austria_changes,
    aes(
      fill = bult_area_change
    ),    color = "white",
    size = .15
  ) +
  scale_fill_gradient2(
    name = "Built Area Change in Austria comparing 2018 and 2023",
    midpoint = 0,
    mid = "yellow",
    high = "red",
    low = "green"
  ) +
  guides(
    fill = guide_colorbar(
      direction = "horizontal",
      barheight = unit(1.5, "mm"),
      barwidth = unit(20, "mm"),
      title.position = "top",
      label.position = "bottom",
      title.hjust = .5,
      label.hjust = .5,
      nrow = 1,
      byrow = T
    )
  ) +
  theme_void() +
  theme(
    legend.position = "top",
    plot.margin = unit(
      c(
        t = 0, b = 0,
        r = 0, l = 0
      ), "lines"
    )
  )
built_map_change_austria

Conclusion

This submission mainly aims to show possibilities in regard to working with terra and sf packages. I didn’t expect it would essentially take such long time to come up with these visualizations, but my main problem was that the documentation regarding these packages is very much limited and only source of information I was able to find were the official documentation pages of the packages. As it seems not many people are working with them and getting along with terra is not the easiest task.

I wanted to get to know these packages as I think they are such a powerful tools for visualizing and inspecting GeoTIF files, so I pointed my time and effort into this direction, not analysing the cvs files I mentioned in our previous conversations. I think these visualizations are meaningful and could be expanded in many more ways. Now that I know how to work GeoTIF data it would be also possible to make meaningful maps for whole continent of Europe exploring diverse data sets that can be found on Global Forest Watch. Combining GeoTIF files for the entire globe is still complicated, but as we mentioned in our previous conversation I will give my best to replicate the tree growth map for the whole globe.

As I stated through out this submission, everything done here was for simple visualization purposes and to see if these things are possible at all. When I discovered Land Cover Data Set I thought i would be perfect for something like this, since it contains information over multiple years and factors. In the further work I would also like to compare diverse factors such as human built areas and human built copes with our variable of interest (trees), to see if there is strong negative correlation between them. I also think it would be good to make an animation of tree cover change over the years provided by Land Cover Data Set which are 2018 to 2023. It would be certainly possible to do so for a single country, as in this case Austria or Serbia, but in the next few days I would try doing so for the whole continent of Europe to see if it’s possible (around 40 tiles - here each tile was around 150Mb).

Since we also have cvs data sets with which we will be working on country base, I found a data set with monthly burned areas world wide: Global Monthly Burned Area [2002 - 2023] (download link) which is in the same form as our data sets, so it would be well suited for the task. Besides that Global Forest Watch features similar data over the years which would also make it possible to use it in our case. I still haven’t downloaded the data and taken a proper look at it, but it can be downloaded from here.

On Global Forest Watch there is also Data Set named Tree Cover Loss by Dominant Driver 2023 (download here). It seems that its resolution isn’t the crispest, but I will take a proper look at it to see if it’s usable in our case.

I also found data about Global Mining, which contains polygons of activities related to mining in polygonal and GeoTIF format, which can be meaningful for our visualizations. It would also be possible to check the correlation between mining activities and deforestation. Data can be found here.